# Run regressions for foot traffic analysis


# Setup ----
options(scipen = 999)
library(caret)
library(tidyverse)
library(lfe)
library(stargazer)
library(cragg)
library(xtable)
library(broom)
library(lubridate)
library(RColorBrewer)



load("/02_Data/FL_visitor_1804_1910_tract.RData") # NOTE: Proprietary data, values of all variables have been replaced with * 
load("/02_Data/tract_flood_FL.RData")             # NOTE: Proprietary data, values of all variables have been replaced with * 
load("/02_Data/NFIP_Michael_Claim.RData")
load("/02_Data/tract2017.RData")
load("/02_Data/NFIPTakeup.RData")
load("/02_Data/FloodZone.RData")
load("/02_Data/SavingsAccount.RData")
load("/02_Data/SCI_v3.RData")


# tract-level data as of 2017 --------------------------------------------------
tract_data = tract2017 %>%
  left_join(tract_flood_FL, by=c("GEOID"="FIPS_tract")) %>%
  left_join(NFIPTakeup, by="GEOID") %>%
  left_join(FloodZone, by="GEOID") %>%
  left_join(SavingsAccount, by="GEOID") %>%
  left_join(SCI_v3, by=c("GEOID"="tract")) %>%
  mutate(FloodDepth_ft_tract_dummy = ifelse(FloodDepth_ft_tract>0, 1, 0))



# tract-business-time level data, April 2018 to October 2019 -------------------
FL_visitor_1804_1910 = FL_visitor_1804_1910_tract %>%
  left_join(NFIP_Michael_Claim, by=c("visitor_home_tract"="CensusTract", "month")) %>%
  left_join(tract_data, by=c("visitor_home_tract"="GEOID")) %>%
  mutate(
    FloodDepth_ft_tract_dummy = ifelse(FloodDepth_ft_tract>0, 1, 0),
    ClaimAmount               = ifelse(is.na(ClaimAmount), 0, ClaimAmount),
    WAvgPaidWeek              = ifelse(is.na(WAvgPaidWeek), 0, WAvgPaidWeek),
    month                     = factor(month)
  )

date.binarized = predict(dummyVars(formula = "~ month", FL_visitor_1804_1910, sep = '_'), FL_visitor_1804_1910)        # the "dummyVars" formula in the "caret" package can be used to binarize the categotical variable automatically
colnames(date.binarized) = substr(colnames(date.binarized)[1:length(colnames(date.binarized))], 7, nchar(colnames(date.binarized)[1:length(colnames(date.binarized))]))

FL_visitor_1804_1910 = cbind.data.frame(FL_visitor_1804_1910, date.binarized)
rm(date.binarized)





# ---------------------------------------
# ------ Propensity score matching ------
# ---------------------------------------
tract_data$FloodInsTakeup_High = ifelse(tract_data$FloodInsTakeup>0.5, 1, 0)

library(MatchIt)

mod_match <- matchit(FloodInsTakeup_High ~ log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                       ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                       ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract, 
                     method = "nearest", 
                     distance = "glm",
                     data = tract_data[tract_data$FloodDepth_ft_tract>0,])

tract_FL_flooded_match = match.data(mod_match)



# Matched sample --------------------------------------------------------------
FL_visitor_1804_1910_matched = FL_visitor_1804_1910 %>%                 
  filter(visitor_home_tract %in% tract_FL_flooded_match$GEOID)











################################################################################
# Figure 2: Event Study for Baseline
################################################################################
a1_Event = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + FloodDepth_ft_tract:FloodInsTakeup + 
                                                     
                                                   FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                                   log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                                   ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                                   ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):(`18-04`+`18-05`+`18-06`+`18-07`+`18-08`+
                                                                                                             `18-10`+`18-11`+`18-12`+`19-01`+`19-02`+
                                                                                                             `19-03`+`19-04`+`19-05`+`19-06`+`19-07`+
                                                                                                             `19-08`+`19-09`+`19-10`) | 
            factor(visitor_home_tract_placekey) + factor(month) | 0 | visitor_home_tract, data = FL_visitor_1804_1910)









# Extracting coefficients and confidence intervals from the regressions
i1_tidy = tidy(a1_Event, conf.int = TRUE)

print(i1_tidy, n = 50) # Values of interest are in rows 39:49

# Only keeping values of interest - coef. of logged flood depth
i1_a <- i1_tidy[1:5, -1]
i1_a[nrow(i1_a)+1,] = 0
i1_b <- i1_tidy[6:18, -1]
i1_vals_temp <- bind_rows(i1_a, i1_b)





# Only keeping values of interest - coef. of logged flood depth
i2_a <- i1_tidy[541:545, -1]
i2_a[nrow(i2_a)+1,] = 0
i2_b <- i1_tidy[546:558, -1]
i2_vals_temp <- bind_rows(i2_a, i2_b)





# Creating date sequence for axis
start_2018 <- as.Date("2018/04/01")
end_2019 <- as.Date("2019/10/01")
Dates_label <- seq.Date(start_2018, end_2019, by = "month")

Month_label <- format(Dates_label, "%Y-%m")



# Figure 3: Panel C Balance (Non-Borrowers vs. Borrowers) --------------
x_seq <- seq(1, 19, 1) # Determines how many x axis labels will be included


# Monthly Balances: Non-borrowers -------------------------------------
i1_vals_temp %>% 
  mutate(Time = c(1:nrow(i1_vals_temp))) %>%
  ggplot(aes(Time, estimate))+
  geom_segment(aes(x = 0, y = 0,  xend = 20, yend = 0), color = "darkgray")+
  geom_vline(xintercept = 7, alpha = 0.5, linetype = "dashed")+ # Lines marking gap in data
  geom_point( size = 5, alpha = 0.8)+ # Point estimates
  geom_line(linewidth = 1, alpha = 0.8) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2)+
  xlab("Event Month (relative to Hurricane Michael)")+
  ylab("Coeff. EventMonth * FloodDepth")+
  coord_cartesian(x = c(1,19), y = c(-2.5, 4))+
  theme_classic()+
  scale_x_continuous(labels = seq(-6,12), breaks = x_seq)+ # Adding month values to axis
  #scale_x_continuous(labels = Month_label[x_seq], breaks = x_seq)+ # Adding month values to axis
  theme(legend.position = "none",
        text = element_text(size = 16)
  ) 

ggsave("/03_Results/01_Figures/Figure2A.jpeg", width = 10, height = 6, dpi=600)


# Monthly Balances: Non-borrowers -------------------------------------
i2_vals_temp %>% 
  mutate(Time = c(1:nrow(i2_vals_temp))) %>%
  ggplot(aes(Time, estimate))+
  geom_segment(aes(x = 0, y = 0,  xend = 20, yend = 0), color = "darkgray")+
  geom_vline(xintercept = 7, alpha = 0.5, linetype = "dashed")+ # Lines marking gap in data
  geom_point(size = 5, alpha = 0.8)+ # Point estimates
  geom_line(linewidth = 1, alpha = 0.8) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.1, fill = "red")+
  xlab("Event Month (relative to Hurricane Michael)")+
  ylab("Coeff. EventMonth * FloodDepth * FloodInsTakeup")+
  coord_cartesian(x = c(1,19), y = c(-2.5, 4))+
  theme_classic()+
  scale_x_continuous(labels = seq(-6,12), breaks = x_seq)+ # Adding month values to axis
  #scale_x_continuous(labels = Month_label[x_seq], breaks = x_seq)+ # Adding month values to axis
  theme(legend.position = "none",
        text = element_text(size = 16)
  ) 


ggsave("/03_Results/01_Figures/Figure2B.jpeg", width = 10, height = 6, dpi=600)





################################################################################
# Figure 3: Event Study for Matched Sample & IV Regression
################################################################################
FL_visitor_1804_1910 = FL_visitor_1804_1910 %>% 
  mutate(
    FloodDepth_FloodInsPtg_1804 = FloodDepth_ft_tract*FloodInsTakeup*`18-04`,
    FloodDepth_FloodInsPtg_1805 = FloodDepth_ft_tract*FloodInsTakeup*`18-05`,
    FloodDepth_FloodInsPtg_1806 = FloodDepth_ft_tract*FloodInsTakeup*`18-06`,
    FloodDepth_FloodInsPtg_1807 = FloodDepth_ft_tract*FloodInsTakeup*`18-07`,
    FloodDepth_FloodInsPtg_1808 = FloodDepth_ft_tract*FloodInsTakeup*`18-08`,
    FloodDepth_FloodInsPtg_1810 = FloodDepth_ft_tract*FloodInsTakeup*`18-10`,
    FloodDepth_FloodInsPtg_1811 = FloodDepth_ft_tract*FloodInsTakeup*`18-11`,
    FloodDepth_FloodInsPtg_1812 = FloodDepth_ft_tract*FloodInsTakeup*`18-12`,
    FloodDepth_FloodInsPtg_1901 = FloodDepth_ft_tract*FloodInsTakeup*`19-01`,
    FloodDepth_FloodInsPtg_1902 = FloodDepth_ft_tract*FloodInsTakeup*`19-02`,
    FloodDepth_FloodInsPtg_1903 = FloodDepth_ft_tract*FloodInsTakeup*`19-03`,
    FloodDepth_FloodInsPtg_1904 = FloodDepth_ft_tract*FloodInsTakeup*`19-04`,
    FloodDepth_FloodInsPtg_1905 = FloodDepth_ft_tract*FloodInsTakeup*`19-05`,
    FloodDepth_FloodInsPtg_1906 = FloodDepth_ft_tract*FloodInsTakeup*`19-06`,
    FloodDepth_FloodInsPtg_1907 = FloodDepth_ft_tract*FloodInsTakeup*`19-07`,
    FloodDepth_FloodInsPtg_1908 = FloodDepth_ft_tract*FloodInsTakeup*`19-08`,
    FloodDepth_FloodInsPtg_1909 = FloodDepth_ft_tract*FloodInsTakeup*`19-09`,
    FloodDepth_FloodInsPtg_1910 = FloodDepth_ft_tract*FloodInsTakeup*`19-10`,
    
    FloodDepth_FloodInsPtg_IV_1804 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-04`,
    FloodDepth_FloodInsPtg_IV_1805 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-05`,
    FloodDepth_FloodInsPtg_IV_1806 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-06`,
    FloodDepth_FloodInsPtg_IV_1807 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-07`,
    FloodDepth_FloodInsPtg_IV_1808 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-08`,
    FloodDepth_FloodInsPtg_IV_1810 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-10`,
    FloodDepth_FloodInsPtg_IV_1811 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-11`,
    FloodDepth_FloodInsPtg_IV_1812 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`18-12`,
    FloodDepth_FloodInsPtg_IV_1901 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-01`,
    FloodDepth_FloodInsPtg_IV_1902 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-02`,
    FloodDepth_FloodInsPtg_IV_1903 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-03`,
    FloodDepth_FloodInsPtg_IV_1904 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-04`,
    FloodDepth_FloodInsPtg_IV_1905 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-05`,
    FloodDepth_FloodInsPtg_IV_1906 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-06`,
    FloodDepth_FloodInsPtg_IV_1907 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-07`,
    FloodDepth_FloodInsPtg_IV_1908 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-08`,
    FloodDepth_FloodInsPtg_IV_1909 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-09`,
    FloodDepth_FloodInsPtg_IV_1910 = FloodDepth_ft_tract*log(SCI_Flooded_Distant)*`19-10`
  )








a3_Event = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + FloodDepth_ft_tract:I(FloodInsTakeup>0.5) + 
                                               
                                                     FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                                     log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                                     ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                                     ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):(`18-04`+`18-05`+`18-06`+`18-07`+`18-08`+
                                                                                                           `18-10`+`18-11`+`18-12`+`19-01`+`19-02`+
                                                                                                           `19-03`+`19-04`+`19-05`+`19-06`+`19-07`+
                                                                                                           `19-08`+`19-09`+`19-10`) | 
            
            factor(visitor_home_tract_placekey) + factor(month) | 0 | visitor_home_tract, data = FL_visitor_1804_1910_matched)

a4_Event = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + 
                                               
                                               FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                               log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                               ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                               ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):(`18-04`+`18-05`+`18-06`+`18-07`+`18-08`+
                                                                                                           `18-10`+`18-11`+`18-12`+`19-01`+`19-02`+
                                                                                                           `19-03`+`19-04`+`19-05`+`19-06`+`19-07`+
                                                                                                           `19-08`+`19-09`+`19-10`) | 
            
            factor(visitor_home_tract_placekey) + factor(month) | 
            (FloodDepth_FloodInsPtg_1804|FloodDepth_FloodInsPtg_1805|FloodDepth_FloodInsPtg_1806|FloodDepth_FloodInsPtg_1807|FloodDepth_FloodInsPtg_1808|FloodDepth_FloodInsPtg_1810|FloodDepth_FloodInsPtg_1811|FloodDepth_FloodInsPtg_1812|FloodDepth_FloodInsPtg_1901|FloodDepth_FloodInsPtg_1902|FloodDepth_FloodInsPtg_1903|FloodDepth_FloodInsPtg_1904|FloodDepth_FloodInsPtg_1905|FloodDepth_FloodInsPtg_1906|FloodDepth_FloodInsPtg_1907|FloodDepth_FloodInsPtg_1908|FloodDepth_FloodInsPtg_1909|FloodDepth_FloodInsPtg_1910 ~ 
                FloodDepth_FloodInsPtg_IV_1804+FloodDepth_FloodInsPtg_IV_1805+FloodDepth_FloodInsPtg_IV_1806+FloodDepth_FloodInsPtg_IV_1807+FloodDepth_FloodInsPtg_IV_1808+FloodDepth_FloodInsPtg_IV_1810+FloodDepth_FloodInsPtg_IV_1811+FloodDepth_FloodInsPtg_IV_1812+FloodDepth_FloodInsPtg_IV_1901+FloodDepth_FloodInsPtg_IV_1902+FloodDepth_FloodInsPtg_IV_1903+FloodDepth_FloodInsPtg_IV_1904+FloodDepth_FloodInsPtg_IV_1905+FloodDepth_FloodInsPtg_IV_1906+FloodDepth_FloodInsPtg_IV_1907+FloodDepth_FloodInsPtg_IV_1908+FloodDepth_FloodInsPtg_IV_1909+FloodDepth_FloodInsPtg_IV_1910) | 
            visitor_home_tract, data = FL_visitor_1804_1910)










# Extracting coefficients and confidence intervals from the regressions
i3_tidy = tidy(a3_Event, conf.int = TRUE)

print(i3_tidy, n = 540) # Values of interest are in rows 39:49

# Only keeping values of interest - coef. of logged flood depth
i3_a <- i3_tidy[523:527, -1]
i3_a[nrow(i3_a)+1,] = 0
i3_b <- i3_tidy[528:540, -1]
i3_vals_temp <- bind_rows(i3_a, i3_b)



i4_tidy = tidy(a4_Event, conf.int = TRUE)

print(i4_tidy, n = 540) # Values of interest are in rows 39:49


# Only keeping values of interest - coef. of logged flood depth
i4_a <- i4_tidy[541:545, -1]
i4_a[nrow(i4_a)+1,] = 0
i4_b <- i4_tidy[546:558, -1]
i4_vals_temp <- bind_rows(i4_a, i4_b)





# Creating date sequence for axis
start_2018 <- as.Date("2018/04/01")
end_2019 <- as.Date("2019/10/01")
Dates_label <- seq.Date(start_2018, end_2019, by = "month")

Month_label <- format(Dates_label, "%Y-%m")



# Figure 3: Panel C Balance (Non-Borrowers vs. Borrowers) --------------
x_seq <- seq(1, 19, 1) # Determines how many x axis labels will be included



# Monthly Balances: Non-borrowers -------------------------------------
i3_vals_temp %>% 
  mutate(Time = c(1:nrow(i3_vals_temp))) %>%
  ggplot(aes(Time, estimate))+
  geom_segment(aes(x = 0, y = 0,  xend = 20, yend = 0), color = "darkgray")+
  geom_vline(xintercept = 7, alpha = 0.5, linetype = "dashed")+ # Lines marking gap in data
  geom_point( size = 5, alpha = 0.8)+ # Point estimates
  geom_line(linewidth = 1, alpha = 0.8) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.1, fill = "red")+
  xlab("Event Month (relative to Hurricane Michael)")+
  ylab("Coeff. EventMonth * FloodDepth * HighFloodInsTakeup")+
  coord_cartesian(x = c(1,19), y = c(-0.7, 2.6))+
  theme_classic()+
  scale_x_continuous(labels = seq(-6,12), breaks = x_seq)+ # Adding month values to axis
  #scale_x_continuous(labels = Month_label[x_seq], breaks = x_seq)+ # Adding month values to axis
  theme(legend.position = "none",
        text = element_text(size = 16)
  ) 

ggsave("/03_Results/01_Figures/Figure3A.jpeg", width = 10, height = 6, dpi=600)


# Monthly Balances: Non-borrowers -------------------------------------
i4_vals_temp %>% 
  mutate(Time = c(1:nrow(i2_vals_temp))) %>%
  ggplot(aes(Time, estimate))+
  geom_segment(aes(x = 0, y = 0,  xend = 20, yend = 0), color = "darkgray")+
  geom_vline(xintercept = 7, alpha = 0.5, linetype = "dashed")+ # Lines marking gap in data
  geom_point(size = 5, alpha = 0.8)+ # Point estimates
  geom_line(linewidth = 1, alpha = 0.8) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.1, fill = "red")+
  xlab("Event Month (relative to Hurricane Michael)")+
  ylab("Coeff. EventMonth * FloodDepth * FloodInsTakeup")+
  coord_cartesian(x = c(1,19), y = c(-2.5, 10))+
  theme_classic()+
  scale_x_continuous(labels = seq(-6,12), breaks = x_seq)+ # Adding month values to axis
  #scale_x_continuous(labels = Month_label[x_seq], breaks = x_seq)+ # Adding month values to axis
  theme(legend.position = "none",
        text = element_text(size = 16)
  ) 


ggsave("/03_Results/01_Figures/Figure3B.jpeg", width = 10, height = 6, dpi=600)



